Introduction à l'Inférence Bayesienne avec Python et PyMC3

Club de lecture Octobre 2018

Alexis Fortin-Côté

How often have I said to you that when you have eliminated the impossible, whatever remains, however improbable, must be the truth?

-Sherlock Holmes, The Sign of the Four Chap. 6, p. 111, (1890)

A partir d'une présentation de @AustinRochford

original notebook

Inférence Bayésienne et Analyses fréquentistes

Bayésienne Fréquentistes

Problème avec les p-values

  • Gelman, A., & Geurts, H. M. (2017). The statistical crisis in science: How is it relevant to clinical neuropsychology?. The Clinical Neuropsychologist, 31(6-7), 1000-1014.
  • Baker, M. (2016). 1,500 scientists lift the lid on reproducibility. Nature News, 533(7604), 452.
  • Gelman, A. (2016). The problems with p-values are not just with p-values. The American Statistician, 70.
  • Button, K. S., Ioannidis, J. P., Mokrysz, C., Nosek, B. A., Flint, J., Robinson, E. S., & Munafò, M. R. (2013). Power failure: why small sample size undermines the reliability of neuroscience. Nature - Reviews Neuroscience, 14(5), 365.
  • Simmons, J., Nelson, L., & Simonsohn, U. (2011). False-positive psychology: Undisclosed fexibility in data collection and analysis allow presenting anything as signifcant. Psychological Science, 22, 1359–1366.
  • Schatz, P., Jay, K. A., McComb, J., & McLaughlin, J. R. (2005). Misuse of statistical tests in archives of clinical neuropsychology publications. Archives of Clinical Neuropsychology, 20, 1053–1059
  • Gelman, A. (2005). Analysis of variance—why it is more important than ever. The annals of statistics, 33(1), 1-53.
  • Bezeau, S., & Graves, R. (2001). Statistical power and effect sizes of clinical neuropsychology research. Journal of Clinical and Experimental Neuropsychology, 23(3), 399-406.
  • Meehl, P. E. (1990). Why summaries of research on psychological theories are often uninterpretable. Psychological Reports, 66, 195–244

Inférence Bayesienne

Exemple classique:

Supposons un test sanguin permettant de détecter le vampirimse avec un taux de succès de 95% Pr(positive|vampire) = 0.95. C'est donc un test précis, mais qui implique son lot de faux-positif. Un pourcent du temps le test sera positif pour une personne normale, donc Pr(positive|mortel)=0.01. La dernière donnée connue est qu'une personne sur mille est un vampire dans la population générale, donc Pr(vampire)=0.001

Probabilité d'être un vampire si on test positif? 95%? 1%?

Probabilité conditionnelle

En théorie des probabilités, une probabilité conditionnelle est la probabilité d'un événement sachant qu'un autre événement a eu lieu.

$$ \begin{align*} P(A\ |\ B) & = \textrm{La probabilité que } A \textrm{ arrive si on sait que } B \textrm{ est arrivé} \\ & = \frac{P(A \textrm{ et } B)}{P(B)}. \end{align*} $$

Notre question,

Quelle est la probabilité qu'une personne soit un vampire si elle est testée positive?

$$ \begin{align*} P(vampire) & = 0.001 \\ P(\textrm{Test } +\ |\ vampire) & = 0.95 \\ P(\textrm{Test } +\ |\ mortel) & = 0.01 \\ \\ P(vampire\ |\ \textrm{Test } +) & =\ \textbf{?} \end{align*} $$

Théorem de Bayes

Le théorem de Bayes nous montre comment passer de $P(B\ |\ A)$ à $P(A\ |\ B)$.

<img src='https://upload.wikimedia.org/wikipedia/commons/d/d4/Thomas_Bayes.gif' width=300> $$\large P(A\ |\ B) = \frac{P(B\ |\ A)\ P(A)}{P(B)}$$

Le terme P(A) est la probabilité a priori de A. Elle est « antérieure » au sens qu’elle précède toute information sur B. Le terme P(A|B) est appelée la probabilité conditionnelle de A sachant B . Elle est « postérieure », au sens qu’elle dépend directement de B. Le terme P(B|A), pour un B connu, est appelé la fonction de vraisemblance de A. De même, le terme P(B) est appelé la probabilité marginale ou a priori de B.

$$\color{red}{P(A\ |\ B)} = \frac{\color{blue}{P(B\ |\ A)}\ \color{orange}{P(A)}}{\color{green}{P(B)}}.$$

Pour plusieurs modèles la la probabilité marginale est impossible ou difficilement calculable analytiquement.

Pour le cas qui nous intéresse, on peut calculer analytiquement la probabilité avec le théorème de Bayes et la Formule des probabilités totales

$$Pr(vampire|+)=\frac{Pr(+|vampire)Pr(vampire)}{Pr(+)}$$
$$ Pr(+)=Pr(+|vampire)Pr(vampire)+Pr(+|mortel)(1−Pr(vampire)) $$
$$ \begin{align*} P(vampire) & = 0.001 \\ P(+\ |\ vampire) & = 0.95 \\ P(+\ |\ mortel) & = 0.01 \\ \\ P(vampire\ |\ \textrm{Test } +) & = \frac{0.95*0.001}{0.95*0.001 + 0.01*(1-0.001)} \end{align*} $$
In [1]:
(0.95*0.001)/(0.95*0.001 + 0.01 * (1-0.001))
Out[1]:
0.08683729433272395

C'est un exemple classique, mais en soi pas particuilièrement Bayésien

Il existe une façon plus intuitive de poser problème

  1. Dans une population de 100,000 personnes, 100 sont vampires.

  2. De ces 100 vampires, 95 vont tester positif au test.

  3. Des 99,900 mortels, 999 vont tester positif au test.

Donc si on teste les 100,000 personnes, quelle est la proportion de ceux qui teste positif d'être des vampires? La question posée ainsi semble plus facile à répondre.

On a juste à compter le nombre de test positifs: 95 + 999 = 1094. On sait que dans le lot 95 sont des vampires:

Pr(vampire|+)=95/1094≈0.087

Programmation probabilistique pour inférence bayesienne

L'exemple précédant était déjà complexe à exprimer et pouvait se résoudre analytiquement parce qu'il était très simple.

Résolvons le même problème avec Python

On connait la prévalence de la maladie dans la population générale, soit 1/1000

La distribution de Bernoulli donne la probabilité binaire (pile ou face biaisé). Si $X \sim \textrm{Bernoulli}(p),$

$$ \begin{align*} P(X = 1) & = p \\ P(X = 0) & = 1 - p. \end{align*} $$
In [10]:
import pymc3 as pm

with pm.Model() as disease_model:
    has_disease = pm.Bernoulli('est_vampire', 0.001)

Si la personne est vampire, il y a 95% de chance de tester positif. Sinon, il y a quand même 1% de chance de tester positif

In [11]:
with disease_model:
    p_test_pos = has_disease * 0.95 + (1 - has_disease) * 0.01

La personne dans l'exemple a testé positive alors:

In [12]:
with disease_model:
    test_pos = pm.Bernoulli('test_pos', p_test_pos, observed=1)

Quelle est la probabilité que cette personne ait la maladie

In [13]:
with disease_model:
    disease_trace = pm.sample(draws=10000, random_seed=SEED, cores=1)
In [14]:
disease_trace['est_vampire']
Out[14]:
array([0, 0, 0, ..., 0, 0, 1], dtype=int64)
In [15]:
disease_trace['est_vampire'].mean()
Out[15]:
0.0862

Ce qu'on a calculé tantôt:

In [16]:
(0.95*0.001)/(0.95*0.001 + 0.01 * (1-0.001))
Out[16]:
0.08683729433272395

Méthodes de Monte Carlo

Les méthodes de Monte Carlo utilisent un échantillionnage aléatoire pour aproximer des quantités difficile a obtenir analytiquement. Ils sont un outil essentiel pour l'inférence Bayesienne où la distribution marginale est habituellement difficile à obtenir.

Un exemple

On génère 5000 points distribués aléatoirement dans un carré unitaire

In [17]:
N = 5000

x, y = np.random.uniform(0, 1, size=(2, N))
In [19]:
fig
Out[19]:

En comptant le nombre de point qui tombent dans le quatier de cercle de rayon 1, on obtient une approximation de l'aire dans de quartier de cercle, qui est en réalité de $\frac{\pi}{4}$.

In [20]:
in_circle = x**2 + y**2 <= 1
In [22]:
fig
Out[22]:
In [23]:
4 * in_circle.mean()
Out[23]:
3.1672

Historique des méthodes de Monte Carlo

Les méthodes de Monte Carlo ont été largement utilisées pour le Projet Manhattan. Sur la photo ci-dessus, on voit des scientifiques du projet Manhattan Stanislaw_Ulam,Richard Feynman et John von Neumann. En travaillant sur le projet Manhattan, Ulam a donné l'une des premières descriptions des algorithmes des Markov Chain Monte Carlo. La librairie de programmation probabiliste Stan est nommé en son honneur.

Étude de cas:

La pièce de monnaie truquée

Nous avons en notre possession une pièce de monnaie. On veut savoir si elle est truquée, c.-à-d. si la probabilité d'obtenir pile ou face est égale

Bayésien = Mise à jour de nos apriori en fonction de nos observations

On peut représenter un lancer de la pièce par un essai de Bernouilli, qui prend la valeur 1 avec la probabilité p.

À priori on ne connait pas la valeur de p. Donc la plausibilité de p est uniforme entre 0 et 1.

In [10]:
x = np.linspace(0, 1, 200)
pdf = stats.beta.pdf(x, 1, 1)
plt.plot(x, pdf)
plt.xlabel('p', fontsize=12)
plt.ylabel('pdf', fontsize=12)
plt.ylim(0, 4.5)
plt.show()

On lance la pièce de monnaie, on obtient face

In [13]:
x = np.linspace(0, 1, 200)
pdf = stats.beta.pdf(x, 2, 1)
plt.plot(x, pdf, label=r'face = {}, pile = {}'.format(1, 0))
plt.legend()
plt.xlabel('p', fontsize=12)
plt.ylabel('pdf', fontsize=12)
plt.ylim(0, 4.5)
plt.show()

On lance la pièce de monnaie, on obtient pile

In [14]:
x = np.linspace(0, 1, 200)
pdf = stats.beta.pdf(x, 2, 2)
plt.plot(x, pdf, label=r'face = {}, pile = {}'.format(1, 1))
plt.legend()
plt.xlabel('p', fontsize=12)
plt.ylabel('pdf', fontsize=12)
plt.ylim(0, 4.5)
plt.show()

On lance la pièce de monnaie, on obtient pile

In [15]:
x = np.linspace(0, 1, 200)
pdf = stats.beta.pdf(x, 2, 3)
plt.plot(x, pdf, label=r'face = {}, pile = {}'.format(1, 2))
plt.legend()
plt.xlabel('p', fontsize=12)
plt.ylabel('pdf', fontsize=12)
plt.ylim(0, 4.5)
plt.show()
In [83]:
HTML(anim.to_html5_video())
Out[83]:

Étude de cas : Le problème de Monty Hall

Le problème de Monty Hall est un célèbre puzzle de probabilité, basé sur l'émission de jeux des années 1960 Let's Make a Deal et nommé d'après son hôte d'origine. Dans le jeu, un concurrent se voyait présenter trois portes, dont deux contenaient un article de peu ou pas de valeur (par exemple, une chèvre) et l'une d'entre elles contenait un article de très grande valeur (par exemple, une voiture de luxe). Le concurrent devait d'abord deviner quelle porte contenait la voiture de sport. Après la devinette initiale du concurrent, Monty ouvrait l'une des deux autres portes, révélant une chèvre. Monty offrait alors au concurrent la chance de changer leur choix de porte. Le problème de Monty Hall pose la question suivante : le candidat doit-il conserver son choix initial de porte ou le changer ?

Au départ, nous n'avons aucune information sur la porte derrière laquelle se trouve le prix.

In [24]:
with pm.Model() as monty_model:
    prize = pm.DiscreteUniform('prize', 0, 2)

Si nous choisisson la porte 1:

Monty peut ouvrir
Prix derrière Porte 1 Porte 2 Porte 3
Porte 1 Non Oui Oui
Porte 2 Non Non Oui
Porte 3 Non Oui Non
In [25]:
from theano import tensor as tt

with monty_model:
    # Probability that Monty open each door
    p_open = pm.Deterministic('p_open',
                              tt.switch(tt.eq(prize, 0),
                                        # it is behind the first door
                                        np.array([0., 0.5, 0.5]),
                              tt.switch(tt.eq(prize, 1),
                                        # it is behind the second door
                                        np.array([0., 0., 1.]),
                                        # it is behind the third door
                                        np.array([0., 1., 0.]))))

Monty ouvre la porte 3 contenant une chèvre

In [26]:
with monty_model:
    opened = pm.Categorical('opened', p_open, observed=2)

Doit-on changer de porte?

In [27]:
with monty_model:
    monty_trace = pm.sample(draws=10000, random_seed=SEED)
    
monty_df = pm.trace_to_dataframe(monty_trace)
In [29]:
fig
Out[29]:

Donc, oui on devrait changer de porte

En option

We can also resolve the Monty Hall problem with pen and paper, as follows.

Throughout this calculation, all probabilities assume that we have initially chosen the first door. By Bayes' Theorem, the probability that the sportscar is behind door one given that Monty opened door three is

$$P(\textrm{Behind door one}\ |\ \textrm{Opened door three}) = \frac{P(\textrm{Opened door three}\ |\ \textrm{Behind door one})P(\textrm{Behind door one})}{P(\textrm{Opened door three})}.$$

The a priori probability that the prize is behind any of the doors is one third. From the table above, $P(\textrm{Opened door three}\ |\ \textrm{Behind door one}) = \frac{1}{2}$. We calculate $P(\textrm{Opened door three})$ using the law of total probability as follows:

$$ \begin{align*} P(\textrm{Opened door three}) & = P(\textrm{Opened door three}\ |\ \textrm{Behind door one})P(\textrm{Behind door one}) \\ & + P(\textrm{Opened door three}\ |\ \textrm{Behind door two})P(\textrm{Behind door two}) \\ & + P(\textrm{Opened door three}\ |\ \textrm{Behind door three})P(\textrm{Behind door three}) \\ & = \frac{1}{2} \cdot \frac{1}{3} + 1 \cdot \frac{1}{3} + 0 \cdot \frac{1}{3} \\ & = \frac{1}{2}. \end{align*} $$

Therefore

$$ \begin{align*} P(\textrm{Behind door one}\ |\ \textrm{Opened door three}) & = \frac{P(\textrm{Opened door three}\ |\ \textrm{Behind door one})P(\textrm{Behind door one})}{P(\textrm{Opened door three})} \\ & = \frac{\frac{1}{2} \cdot \frac{1}{3}}{\frac{1}{3}} \\ & = \frac{1}{3}. \end{align*}$$

Since $P(\textrm{Behind door three}\ |\ \textrm{Opened door three}) = 0$, because Monty wants the contestant's choice to be suspensful, $P(\textrm{Behind door two}\ |\ \textrm{Opened door three}) = \frac{2}{3}$. Therefore it is correct to switch doors, confirming our computational results.

Étude de cas: Sleep Deprivation

Les données proviennent de la librairie lme4 du language R, qui cite

Gregory Belenky, Nancy J. Wesensten, David R. Thorne, Maria L. Thomas, Helen C. Sing, Daniel P. Redmond, Michael B. Russo and Thomas J. Balkin (2003) Patterns of performance degradation and restoration during sleep restriction and subsequent recovery: a sleep dose-response study. Journal of Sleep Research 12, 1–12.

In [31]:
sleep_df.head()
Out[31]:
reaction days subject reaction_std
0 249.5600 0 308 -0.868968
1 258.7047 1 308 -0.706623
2 250.8006 2 308 -0.846944
3 321.4398 3 308 0.407108
4 356.8519 4 308 1.035777

Dans cette étude, chacun des sujets obtient son temps normal de sommeil lors du premier jour. Leur temps de sommeil est réduit à 3 heures pour les journées subséquentes. La colonnereaction est leur temps de réponse moyen à un nombre de tests. La colonne reaction_std est le résultat de la standardisation (Z-score) reaction sur tous les sujets et tous les jours.

In [33]:
fig
Out[33]:

Chaque sujet devrait avoir un temps de réaction de base qui ne devrait pas trop différer les uns des autres

In [35]:
with pm.Model() as sleep_model:
    μ_α = pm.Normal('μ_α', 0., 5.)
    σ_α = pm.HalfCauchy('σ_α', 5.)
    α = pm.Normal('α', μ_α, σ_α, shape=n_subjects)

Le taux de croissance des temps de réponse en fonction du nombre de jours de privation de sommeil différe selon chacun des sujets.

In [36]:
with sleep_model:
    μ_β = pm.Normal('μ_β', 0., 5.)
    σ_β = pm.HalfCauchy('σ_β', 5.)
    β = pm.Normal('β', μ_β, σ_β, shape=n_subjects)

La combinaison du temps de réponse de base et du taux de croissance correspondent à nos observations.

In [37]:
with sleep_model:
    μ = pm.Deterministic('μ', α[subject_ix] + β[subject_ix] * days)
    σ = pm.HalfCauchy('σ', 5.)
    obs = pm.Normal('obs', μ, σ, observed=reaction_std)

Ce type de model est connu sous le nom de model linéaire hiérarchique (mixe) "hierarchical (or mixed) linear model", parce que il permet à la pente et à la valeur initiale de variée selon le sujet, mais ajoute un terme de régularisation grâce à un apriori partagé par les sujets. Pour plus d'informations lire Gelman and Hill Data Analysis Using Regression and Multilevel/Hierarchical Models.

In [38]:
N_JOBS = 3
JOB_SEEDS = [SEED + i for i in range(N_JOBS)]

with sleep_model:
    sleep_trace = pm.sample(njobs=N_JOBS, random_seed=JOB_SEEDS)

Diagnostic de la convergence

PyMC3 met à disposition une suite d'outils de diagnostic de convergence statistcal convergence diagnostics pour s'assurer que les échantillions soient une bonne approximation de la vrai distribution postérieure.

Les tracés d'énergie et la fraction bayésienne d'informations manquantes (BFMI) sont deux diagnostics pertinents pour l'échantillonnage Monte Carlo Hamiltonien, que nous avons utilisé dans cet exemple. Si les deux distributions dans le tracé d'énergie diffèrent de façon significative (surtout dans les queues), l'échantillonnage n'a pas été très efficace. La fraction bayésienne des informations manquantes quantifie cette différence avec un nombre compris entre zéro et un. Un BFMI proche de 1 est préférable, et un BFMI inférieur à 0,2 indique des problèmes d'efficacité.

Pour plus d'information : Michael Betancourt A Conceptual Introduction to Hamiltonian Monte Carlo et Robust Statistical Workflow with PyStan.

In [39]:
ax = pm.energyplot(sleep_trace, legend=False)
ax.set_title("BFMI = {:.3f}".format(pm.bfmi(sleep_trace)));

Contrairement aux tracés énergétiques et aux BFMI, qui sont spécifiques aux algorithmes de Monte Carlo Hamiltonien, la statistique Gelman-Rubin est applicable à n'importe quel algorithme MCMC, à condition que plusieurs chaînes aient été échantillonnées. Une statistique de Gelman-Rubin près de 1 est préférable, et les valeurs inférieures à 1,1 sont généralement considérées comme une indication de convergence.

In [40]:
max(np.max(gr_values) for gr_values in pm.gelman_rubin(sleep_trace).values())
Out[40]:
1.0042896858571944

Prediction

En inférence Bayesienne les prédictions sont effectuées en échantillionnant la distribution postiérieur, qui est la distribution des observations futures possibles sachant les observations antérieures.

In [41]:
with sleep_model:
    pp_sleep_trace = pm.sample_ppc(sleep_trace)
In [43]:
pp_df.head()
Out[43]:
reaction days subject reaction_std pp_0 pp_1 pp_2 pp_3 pp_4 pp_5 ... pp_490 pp_491 pp_492 pp_493 pp_494 pp_495 pp_496 pp_497 pp_498 pp_499
0 249.5600 0 308 -0.868968 -0.890141 0.514229 -1.814552 -0.418134 -0.797194 -0.459675 ... -1.766223 -1.650690 -1.190322 -1.697140 -0.994287 -0.444406 -0.749693 -0.158537 0.164312 -0.249929
1 258.7047 1 308 -0.706623 -0.725305 -0.233810 -0.445534 -0.191114 0.007865 -1.152587 ... -1.001274 -0.505686 -0.054972 -0.831341 -0.160696 -0.148323 0.359490 -0.306998 -0.083120 -0.496599
2 250.8006 2 308 -0.846944 -0.296069 -0.014863 0.204777 -0.458021 -0.294116 -0.730840 ... -0.834290 -0.482197 0.544785 -0.101860 0.136675 0.654274 0.485466 -0.432618 0.295079 -0.636435
3 321.4398 3 308 0.407108 0.412651 -0.444827 -0.145654 -0.476851 0.622550 0.807614 ... 0.405701 -0.136766 0.364942 0.469502 -0.292710 0.062846 0.470606 0.850621 0.085148 0.335992
4 356.8519 4 308 1.035777 0.594311 0.634889 -0.028154 -0.046365 0.835659 0.409120 ... 0.531687 0.129832 0.884408 0.907722 0.276634 0.351549 1.066166 0.641605 0.340408 -0.356312

5 rows × 504 columns

In [46]:
grid.fig
Out[46]:

Étude de cas: EEVEE

Statitique de l'échantillon

In [12]:
df.groupby('participant').mean().age_participant.plot.hist()
plt.title('age des paticipants')
plt.show()
In [17]:
df.astype(float).groupby('participant').mean().describe().iloc[:, 4:]
Out[17]:
participant_masculin age_participant travail_domaine jeux_video anime
count 121.000000 121.000000 121.000000 121.000000 121.000000
mean 0.280992 29.462810 0.008264 0.876033 0.760331
std 0.451352 10.369058 0.090909 0.330914 0.428657
min 0.000000 19.000000 0.000000 0.000000 0.000000
25% 0.000000 22.000000 0.000000 1.000000 1.000000
50% 0.000000 26.000000 0.000000 1.000000 1.000000
75% 1.000000 33.000000 0.000000 1.000000 1.000000
max 1.000000 70.000000 1.000000 1.000000 1.000000

Statistique des avatars

In [207]:
show_avatar()

Concentrons-nous sur le réalisme

In [11]:
df_realisme = df.dropna(subset=['realisme'])
df_realisme[['realisme', 'avatar']].groupby('avatar').describe()
Out[11]:
realisme
count mean std min 25% 50% 75% max
avatar
1 242.0 51.260331 24.765171 0.0 35.25 53.5 69.00 100.0
2 242.0 49.623967 25.616448 0.0 33.25 53.0 67.00 100.0
3 242.0 74.859504 20.444564 8.0 64.00 77.0 92.00 100.0
4 242.0 62.714876 22.655944 3.0 48.00 63.5 81.75 100.0
5 242.0 45.698347 23.625798 0.0 31.25 43.0 64.00 99.0
6 242.0 44.880165 23.456106 0.0 33.00 42.0 60.00 100.0
7 242.0 46.438017 25.608088 0.0 28.00 43.5 66.75 100.0
8 242.0 61.962810 22.789190 0.0 46.25 62.0 80.00 100.0
9 242.0 54.913223 24.112218 0.0 38.00 57.0 72.75 100.0
10 242.0 61.504132 21.009432 0.0 50.00 64.0 76.75 100.0
11 242.0 68.541322 20.967282 7.0 57.00 70.0 84.75 100.0
12 242.0 73.685950 20.332315 6.0 62.00 77.0 90.00 100.0
13 242.0 67.330579 22.371028 0.0 57.00 70.0 84.00 100.0
14 242.0 48.954545 24.047369 0.0 32.00 47.0 65.00 100.0
15 242.0 71.590909 21.324844 0.0 59.00 75.0 89.00 100.0
16 242.0 71.537190 19.960975 8.0 60.00 74.0 86.75 100.0

Premier model simple

On crée notre model génératif
Supposons que chaque observation du réalisme pour chacun des avatars provient d'une distribution normale, ayant une moyenne et une variance.

In [12]:
with pm.Model() as model_simple_realisme:
    # realisme
    ## priors
    μ_avatars_realisme = pm.Uniform('μ_avatars_realisme', 0, 100, shape=num_avatar)
    σ_avatars_realisme = pm.HalfCauchy('σ_avatars_realisme', 25, shape=num_avatar)
    
    ## obsevations
    obs_realisme = pm.Normal('obs_realisme', 
                        mu=μ_avatars_realisme[df_realisme.avatar.cat.codes], 
                        sd=σ_avatars_realisme[df_realisme.avatar.cat.codes],
                        observed=df_realisme.realisme)
    
display(model_to_graphviz(model_simple_realisme))
%3 cluster16 16 cluster3,872 3,872 μ_avatars_realisme μ_avatars_realisme ~ Uniform obs_realisme obs_realisme ~ Normal μ_avatars_realisme->obs_realisme σ_avatars_realisme σ_avatars_realisme ~ HalfCauchy σ_avatars_realisme->obs_realisme
In [55]:
with model_simple_realisme:
    trace_simple_real = pm.sample(draws=2000, cores=cores, chains=1)
In [56]:
pm.plots.traceplot(trace_simple_real)
plt.show()
In [57]:
pm.plot_posterior(trace_simple_real, varnames=['μ_avatars_realisme'])
plt.show()
In [58]:
pm.plots.forestplot(trace_simple_real, varnames=['μ_avatars_realisme'])
plt.show()
In [290]:
pm.plots.forestplot(trace_simple_real, varnames=['σ_avatars_realisme'])
plt.show()

On peut maintenant critiquer notre modèle

In [60]:
fig
Out[60]:

On remarque que nos distributions ne concordent pas tellement lorsque les observations sont près de 100

On sous-estime donc le réalisme des avatars les plus réalistes

Utilisons une distribution Normale tronquée

entre 0 et 100

In [13]:
with pm.Model() as model_simple_realisme_tronquee:
    # realisme
    ## priors
    μ_avatars_realisme = pm.Uniform('μ_avatars_realisme', 0, 100, shape=num_avatar)
    σ_avatars_realisme = pm.HalfNormal('σ_avatars_realisme', 25, shape=num_avatar)
    
    ## obsevations
    obs_realisme = pm.TruncatedNormal('obs_realisme', 
                        mu=μ_avatars_realisme[df_realisme.avatar.cat.codes], 
                        sd=σ_avatars_realisme[df_realisme.avatar.cat.codes],
                        upper=100,
                        lower=0,
                        observed=df_realisme.realisme)
    
display(model_to_graphviz(model_simple_realisme_tronquee))
%3 cluster16 16 cluster3,872 3,872 μ_avatars_realisme μ_avatars_realisme ~ Uniform obs_realisme obs_realisme ~ TruncatedNormal μ_avatars_realisme->obs_realisme σ_avatars_realisme σ_avatars_realisme ~ HalfNormal σ_avatars_realisme->obs_realisme
In [63]:
pm.plots.traceplot(trace_model_simple_realisme_tronquee)
plt.show()

On a permis aux moyennes d'être près des limites

In [212]:
pm.plots.forestplot(trace_model_simple_realisme_tronquee, varnames=['μ_avatars_realisme'], chain_spacing=0)
plt.show()

On peut maintenant comparer avec l'ancien modèle

In [66]:
fig
Out[66]:

Comparaison avec le Widely Available Information Criterion

In [67]:
model_simple_realisme_tronquee.name = 'Normale Tronquée'
model_simple_realisme.name = 'Normale'
pm.compare({model_simple_realisme_tronquee:trace_model_simple_realisme_tronquee,
           model_simple_realisme:trace_simple_real})
Out[67]:
WAIC pWAIC dWAIC weight SE dSE var_warn
Normale Tronquée 34349.3 31.28 0 1 63.11 0 0
Normale 35187.6 30.32 838.26 0 84.36 41.76 0

On peut à présent analyser nos inférence pour répondre à notre question:

Quels sont les avatars les plus réalistes?

In [284]:
pd.DataFrame(data={'realisme':trace_model_simple_realisme_tronquee['μ_avatars_realisme'].mean(axis=0), 'avatar':np.arange(1,17)}).sort_values('realisme', ascending=False).set_index('avatar').head()
Out[284]:
realisme
avatar
3 96.772936
12 95.535196
15 94.267306
16 89.654039
13 87.179722

Par contre, ces inférences sont affectés par le biais de notre échantillon, e.g. majoritairement des femmes dans la vingtaine.

Aussi, certains participants peuvent avoir un fort biais ayant un grand impact sur nos réponses.

Essayons de régler ces problèmes

Bâtissons un modèle hierarchique

Permettons à chacun des participants d'être biaisés.
On veut aussi estimer quelle est la variance dans le biais du groupe de participants

In [15]:
display(model_to_graphviz(model_herarchical_real))
%3 cluster121 121 cluster16 16 cluster3,872 3,872 σ_group σ_group ~ HalfCauchy μ_participant_bias_realisme μ_participant_bias_realisme ~ Normal σ_group->μ_participant_bias_realisme obs_realisme obs_realisme ~ TruncatedNormal μ_participant_bias_realisme->obs_realisme μ_avatars_realisme μ_avatars_realisme ~ Uniform μ_avatars_realisme->obs_realisme σ_avatars_realisme σ_avatars_realisme ~ HalfCauchy σ_avatars_realisme->obs_realisme
In [231]:
pm.plots.forestplot(trace_model_herarchical_real, varnames=['μ_avatars_realisme'], rhat=False, chain_spacing=0)
plt.show()
In [232]:
pm.plots.forestplot(trace_model_herarchical_real, varnames=['μ_participant_bias_realisme'], rhat=False, chain_spacing=0)
plt.yticks(plt.yticks()[0], ['' if idx%4 else y for (idx, y) in enumerate(plt.yticks()[1]) ])
plt.show()

On peut toujours comparer avec les autres modèles

In [241]:
pm.densityplot([trace_model_simple_realisme_tronquee, trace_model_herarchical_real], models=['simple', 'herarchique'])
plt.show()
In [72]:
model_simple_realisme_tronquee.name = 'Normale Tronquée'
model_simple_realisme.name = 'Normale'
model_herarchical_real.name = 'Herarchique'

pm.compare({model_simple_realisme_tronquee:trace_model_simple_realisme_tronquee,
           model_simple_realisme:trace_simple_real,
           model_herarchical_real:trace_model_herarchical_real})
Out[72]:
WAIC pWAIC dWAIC weight SE dSE var_warn
Herarchique 32357.7 151.26 0 0.97 91.74 0 1
Normale Tronquée 34349.3 31.28 1991.65 0.03 63.11 89.89 0
Normale 35187.6 30.32 2829.91 0 84.36 108.52 0

Modélisons maintenant en incluant plusieurs prédicteurs

In [17]:
display(graph)
%3 cluster16 16 cluster3,872 3,872 cluster121 121 μ_avatars_realisme μ_avatars_realisme ~ Uniform μ_avatar_bias_realisme μ_avatar_bias_realisme ~ Deterministic μ_avatars_realisme->μ_avatar_bias_realisme σ_avatars_realisme σ_avatars_realisme ~ HalfCauchy obs_realisme obs_realisme ~ TruncatedNormal σ_avatars_realisme->obs_realisme β_avatar_age β_avatar_age ~ Normal β_avatar_age->μ_avatar_bias_realisme β_videogame_par β_videogame_par ~ Normal μ_participant_bias_realisme μ_participant_bias_realisme ~ Deterministic β_videogame_par->μ_participant_bias_realisme β_age_par β_age_par ~ Normal β_age_par->μ_participant_bias_realisme γ_avatar_par_same_genre γ_avatar_par_same_genre ~ Normal γ_avatar_par_same_genre->obs_realisme σ_group σ_group ~ HalfCauchy μ_participant_intercept μ_participant_intercept ~ Normal σ_group->μ_participant_intercept β_travail_par β_travail_par ~ Normal β_travail_par->μ_participant_bias_realisme β_anime_par β_anime_par ~ Normal β_anime_par->μ_participant_bias_realisme β_ismale_participant β_ismale_participant ~ Normal β_ismale_participant->μ_participant_bias_realisme β_avatar_genre β_avatar_genre ~ Normal β_avatar_genre->μ_avatar_bias_realisme μ_participant_intercept->μ_participant_bias_realisme

Les prédicteurs liés au participant

In [291]:
pm.plot_posterior(trace_model_herarchical_real_all_bias_par, varnames=['β_age_par', 'β_ismale_participant', 'β_videogame_par', 'β_anime_par', 'β_travail_par',])
plt.show()

Les prédicteurs liés à l'avatar

In [77]:
pm.plot_posterior(trace_model_herarchical_real_all_bias_par, varnames=['γ_avatar_par_same_genre', 'β_avatar_genre', 'β_avatar_age',])
plt.show()
In [79]:
pm.plots.forestplot(trace_model_herarchical_real_all_bias_par, varnames=['μ_avatars_realisme'], rhat=False, chain_spacing=0)
plt.show()

On peut toujours comparer avec les autres modèles

In [81]:
model_simple_realisme_tronquee.name = 'Normale Tronquée'
model_simple_realisme.name = 'Normale'
model_herarchical_real.name = 'Herarchique'
model_herarchical_real_all_bias_par.name = 'Herarchique + covariables'

pm.compare({model_simple_realisme_tronquee:trace_model_simple_realisme_tronquee,
           model_simple_realisme:trace_simple_real,
           model_herarchical_real:trace_model_herarchical_real,
           model_herarchical_real_all_bias_par:trace_model_herarchical_real_all_bias_par,
           })
Out[81]:
WAIC pWAIC dWAIC weight SE dSE var_warn
Herarchique 32357.7 151.26 0 0.8 91.74 0 1
Herarchique + covariables 32358.7 152.05 1.03 0.17 91.28 3.48 1
Normale Tronquée 34349.3 31.28 1991.65 0.03 63.11 89.89 0
Normale 35187.6 30.32 2829.91 0 84.36 108.52 0

Simulation

En plus d'interpréter les données et de mieux comprendre les interactions entre le modèle et les différentes vairable, on peut utiliser notre modèle pour générer des données et simuler un nouvel échantillion, tout en gardant la notion d'incertitude/plausibilité.

Simulons un échantillon provennant d'une population plus diverse que celle obtenue par les liste de courriels de l'Université

In [25]:
fig
Out[25]:

Finalement on peut comparer quels sont les avatars les plus réalistes en fonction du modèle

In [346]:
pd.concat([df_old.reset_index(), df_new.reset_index()] ,axis=1).head(8).style.hide_index()
Out[346]:
avatar realisme échantillion orig. avatar realisme population variée
3 96.7729 3 88.7247
12 95.5352 15 88.4172
15 94.2673 12 82.2926
16 89.654 16 80.213
13 87.1797 11 77.0147
11 84.2837 13 75.4732
4 72.8789 10 71.8286
8 71.2292 8 69.7687

Pour plus d'info

Statistiques bayésiennes


Plusieurs sont disponible à la bibliothèque ou en ligne avec Ariane 2.0

Échosystème de programmation probabiliste

Probabilistic Programming System Language License Discrete Variable Support Automatic Differentiation/Hamiltonian Monte Carlo Variational Inference
PyMC3 Python Apache V2 Yes Yes Yes
Stan C++, R, Python, ... BSD 3-clause No Yes Yes
Edward Python, ... Apache V2 Yes Yes Yes
BUGS Standalone program, R GPL V2 Yes No No
JAGS Standalone program, R GPL V2 Yes No No
SPSS Standalone program Proprietary ? ? ?

Pas seulement pour les ingénieurs ou les statisticiens!

Plusieurs des développeurs de PyMC3 et Stan proviennent des sciences sociales, de psychologie et neurosciences.